// 
// The model used here is first published in Hu2009_CDC, 
// Simultaneous Regulation of Film Thickness, Surface Roughness
// and Porosity in a Multiscale Thin Film Growth Process
// 
//
#define _CRT_SECURE_NO_WARNINGS
#include "diff_equ_1D_gas_phase_rec.h"
#include "../RateLaw/RateLaw_Kleijn1989.h"

const double VERSION = 0.1;
const double TOLERANCE = 1e-7;

diff_equ_1D_gas_phase_rec::diff_equ_1D_gas_phase_rec(double iK,
								   double iZ_MAX,
								   double iP,
								   double iT_gas,
								   double iT_surface,
								   double iX_BC1,
								   double iX_IC,
								   double it_end,
								   double idz,
								   double idt,
								   double idt_rec,
								   char*  iFileName,
								   int iDetail)
{
	double D;
	double c0;
	double c1;
	int i;
	// Model settings and initial values
	Z_MAX = iZ_MAX;
	K = iK;
	P = iP;
	X_ini = iX_IC;
	Xin_ini = iX_BC1;

	// Simulation settings
	dz = idz;
	t_begin = 0.0;
	t_end = it_end;
	dt_rec = idt_rec;
	dt = idt;
	OutputFileName = iFileName;
	Detail=iDetail;
	if(Detail==1)
	{
		strcpy(FileName2,"");
		strncat(FileName2,OutputFileName,strlen(OutputFileName)-4);
		strcat(FileName2,"_detail.dat");
		fout_detail.open(FileName2);
	}
	// I/O
	x_BC1 = Xin_ini;

	// Internal variables
	znode = int(Z_MAX/dz)+1;
	zARR = new Vec_DP(0.0,znode);
	for(i = 0;i<znode;i++)
	{
		(*zARR)[i] = i*dz;
	}
	t_node = int(t_end/dt_rec)+1;
	t_idx = 0;
	t_arr = new double[t_node];
	for(i=0;i<t_node;i++)
	{
		t_arr[i] = t_begin+i*dt_rec;
	}
	Repeated = 0;
	vecT = new Vec_DP(0.0,znode);
	for(i =0;i<znode-1;i++)
	{
		(*vecT)[i] = iT_gas;
	}
	(*vecT)[znode-1] = iT_surface;
	x = new Vec_DP(0.0,znode-1);
	for(i=0;i<znode-1;i++)
	{
		(*x)[i] = X_ini;
	}
	Xs_sum = new double[t_node];
	Xs2_sum = new double[t_node];
	DepRate_sum = new double[t_node];
	DepRate2_sum = new double[t_node];
	for(i=0;i<t_node;i++)
	{
		Xs_sum[i] = 0.0;
		Xs2_sum[i] = 0.0;
		DepRate_sum[i] = 0.0;
		DepRate2_sum[i] = 0.0;
	}
	// For implcit method
	A = new Mat_DP(0.0,znode-1,znode-1);
	for(i=0;i<znode-1;i++)
	{
		D = diffusivity((*vecT)[i+1],P);
		c0 = -D/(dz*dz);
		c1 = 1/dt+2*D/(dz*dz)+K;

		if(i==0)
		{
			(*A)[i][i] = c1;
			(*A)[i][i+1] = c0;
		}
		else if(i==znode-2)
		{
			(*A)[znode-2][znode-3] = -2*D/(dz*dz);
			(*A)[znode-2][znode-2] = 1/dt+2*D/(dz*dz);
		}
		else
		{
			(*A)[i][i-1] = c0;
			(*A)[i][i]	 = c1;
			(*A)[i][i+1] = c0;
		}
	}
	
	indx = new Vec_INT(0,znode-1);

	LU = new Mat_DP(0.0,znode-1,znode-1);
	(*LU) = (*A);
	ludcmp(*LU,*indx,d);

	b = new Vec_DP(znode-1);

	// Public variables
	time = 0.0;
	Xs = (*x)[znode-2];
	DepRate = ReactionModel(Xs,P,(*vecT)[znode-1]);
}
// ================================================================
diff_equ_1D_gas_phase_rec::~diff_equ_1D_gas_phase_rec()
{
	int i;
	double	Xs_mean;
	double	Xs_var;
	double	Xs_var2;
	double	DepRate_mean;
	double	DepRate_var;
	double	DepRate_var2;

	fout.open(OutputFileName);
	fout << "class: diff_equ_1D_gas_phase_rec v" << VERSION << "\n";
	fout << "Z_max\t= " << Z_MAX << "\n";
	fout << "P\t= " << P << "\n";
	fout << "K\t= " << K << "\n";
	fout << "Ts\t= " << (*vecT)[znode-1] << "\n";
	fout << "X_ini\t= " << X_ini << "\n";
	fout << "Xin_ini\t= " << Xin_ini << "\n";
	fout << "dz\t= " << dz << "\n";
	fout << "t_begin\t= " << t_begin << "\n";
	fout << "t_end = " << t_end << "\n";
	fout << "dt_int = " << dt << "\n";
	fout << "dt_rec = " << dt_rec << "\n";
	fout << "Repeated = " << Repeated << "\n";
	fout << "TOLERANCE = " << TOLERANCE << "\n";
	
	fout << setw(6) << "t(s)" << setw(15) << "<Xs>" << setw(15) << "Var(Xs)" << setw(15) << "<w>(layer/s)" << setw(15) << "Var(w)(layer/s)" << setw(15) << "<Xin>" << setw(15) << "Var(Xin)";
	for(i=0;i<t_node;i++)
	{
		Xs_mean = Xs_sum[i]/Repeated;
		Xs_var2 = (Xs2_sum[i]-Repeated*Xs_mean*Xs_mean)/Repeated;
		if(Xs_var2 < -TOLERANCE)
		{
			cerr << "ERROR@diff_equ_1D_gas_phase_rec::~diff_equ_1D_gas_phase_rec()\n";
			cerr << "i=" << i << " t=" << t_arr[i] << " Xs_var2=" << Xs_var2 << "\n";
			system("pause");
		}
		else if(Xs_var2 >= -TOLERANCE && Xs_var2 < 0.0)
		{
			Xs_var2 = 0.0;
		}
		Xs_var = sqrt(Xs_var2);

		DepRate_mean = DepRate_sum[i]/Repeated;
		DepRate_var2 = (DepRate2_sum[i]-DepRate_mean*DepRate_mean*Repeated)/Repeated;
		if(DepRate_var2 < -TOLERANCE)
		{
			cerr << "ERROR@diff_equ_1D_gas_phase_rec::~diff_equ_1D_gas_phase_rec()\n";
			cerr << "i=" << i << " t=" << t_arr[i] << " DepRate_var2=" << DepRate_var2 << "\n";
			system("pause");
		}
		else if(DepRate_var2 >= -TOLERANCE && DepRate_var2 < 0.0)
		{
			DepRate_var2 = 0.0;
		}
		DepRate_var = sqrt(DepRate_var2);
		
		fout << "\n";
		fout << setw(6) << setprecision(6) << t_arr[i];
		fout << setw(15) << setprecision(6) << Xs_mean;
		fout << setw(15) << setprecision(6) << Xs_var;
		fout << setw(15) << setprecision(6) << DepRate_mean;
		fout << setw(15) << setprecision(6) << DepRate_var;
	}
	fout.close();
	if(Detail==1) fout_detail.close();

	delete zARR;
	delete [] t_arr;
	delete [] Xs_sum;
	delete [] Xs2_sum;
	delete [] DepRate_sum;
	delete [] DepRate2_sum;
	delete vecT;
	delete A;
	delete LU;
	delete indx;
	delete x;
	delete b;
}
// ================================================================
double diff_equ_1D_gas_phase_rec::update(double Xin)
// Use implicit method to solve diffusion equation
//   \frac{\partial X}{\partial t} = D\frac{\partial^2 X}{\partial z^2}-KX
// Discretization
//   \
//   
{
	int i;
	double D2;  // Diffusivity, m^2/s
	double C;   // Concentration, mol/m^3

	if(Xin < 0||Xin>1.0)
	{
		cerr << "ERROR@diff_equ_1D_gas_phase_rec::update():Invalid Xin = " << Xin << "\n";
		system("pause");exit(0);
	}
	x_BC1 = Xin;

	D2 = diffusivity((*vecT)[1],P);

	C  = P/(8.314*(*vecT)[znode-1]);

	(*b)[0] = D2/(dz*dz)*(x_BC1)+1/dt*((*x)[0]);
	for(i=1;i<znode-2;i++)
	{
		(*b)[i] = 1/dt*((*x)[i]);
	}	
	// 2009.11.01 by Xinyu Zhang: Fix the following BIG bug
	// Incorrect one: 
	// (*b)[znode-2] = 1/dt*((*x)[znode-2])-2/dz*(DepRate/5.7423e4);
	(*b)[znode-2] = 1/dt*((*x)[znode-2])-2/dz*(DepRate/5.7423e4/C);

	lubksb(*LU,*indx,*b);
	(*x) = (*b);
	
	Xs = (*x)[znode-2];
	// 2009.07.29 Gangshi: test if Xs is negative but very close to zero
	if ((Xs>-1e-5)&&(Xs<0.0)) Xs=0.0;
	DepRate = ReactionModel(Xs,P,(*vecT)[znode-1]);

	time += dt;
	while(abs((time-t_arr[t_idx])/dt_rec)<0.1)
	{
		Xs_sum[t_idx] += Xs;
		Xs2_sum[t_idx] += Xs*Xs;
		DepRate_sum[t_idx] += DepRate;
		DepRate2_sum[t_idx] += DepRate*DepRate;
		t_idx++;
	}
	return DepRate;
}
// ================================================================
void diff_equ_1D_gas_phase_rec::reset()
{
	Repeated++;
	time = 0.0;
	t_idx = 0;

	for(int i=0;i<znode-1;i++)
	{
		(*x)[i] = X_ini;
	}
	Xs = (*x)[znode-2];

	Xs_sum[t_idx] += Xs;
	Xs2_sum[t_idx] += Xs*Xs;
	DepRate_sum[t_idx] += DepRate;
	DepRate2_sum[t_idx] += DepRate*DepRate;
	t_idx++;
}
// ================================================================
int diff_equ_1D_gas_phase_rec::notfinish()
{
	return (time < t_end) && (t_idx < t_node);
}
// ================================================================
double diff_equ_1D_gas_phase_rec::getXs()
{
	return (*x)[znode-2];
}
// ================================================================
double diff_equ_1D_gas_phase_rec::getXi()
{
	return x_BC1;
}
// ================================================================
double diff_equ_1D_gas_phase_rec::getDepRate()
{
	return DepRate;
}
// ================================================================
void diff_equ_1D_gas_phase_rec::setXi(double iXi)
{
	x_BC1 = iXi;
}

// ================================================================
//
// PRIVATE FUNCTIONS
//
// ================================================================
double diff_equ_1D_gas_phase_rec::diffusivity(double T,double P)
{
	// Unit m^2/s
	return ((-2.90)+(2.06e-2)*T+(2.81e-5)*T*T)/P;
}
// ================================================================

